R Notebook

Load scripts

library(knitr)

source("./Scripts/Settings.R")
source("./Scripts/AusWeather.R")
source("./Scripts/Utility.R")

library(recipes)

aus_map <- function(data) {
    
    map <- with_nameSpace(
        "maps", 
        ggplot2::borders("world", regions = "Australia",  fill = "grey90", colour = "black")
    )
    
    ggplot(data) + map  + xlim(112, 155) + ylim(-44, -10)  # Set y axis limits
}

aus_dark_map <- function(data) {
    
    map <- with_nameSpace(
        "maps", 
        ggplot2::borders("world", regions = "Australia",  fill = "DarkGrey", colour = "black")
    )
    
    ggplot(data) + map  + xlim(112, 155) + ylim(-44, -10)  # Set y axis limits
}

Determine baseline stations

baseline_start <- 1951
baseline_end <- 1980
observation_info <- aus_rainfall_observations_info
minor_period <- aus_monthly_period()
clean_info <- aus_clean_info("RainAnamoly", aus_rainfall_observations_info, yearly_present_target = 350, clean_years_target = 50)

valuation_info <- aus_valuation_info(
    "RainAnamoly", clean_info, baseline_start = baseline_start, baseline_end = baseline_end, 
    minor_period = aus_monthly_period(), 
    scale = TRUE)

baseline_stations <- aus_baseline_stations(valuation_info)
baseline_clean_summaries <- aus_observation_clean_summary(baseline_stations$Station, clean_info) %>% filter(Clean)
baseline_stations

baseline_clean_summaries 

baseline_clean_summaries  %>% 
    unnest(All) %>%
    aus_dark_map() + 
    geom_point(aes(x = Longitude, y = Latitude, color = log2(Mean))) +
    scale_colour_distiller(type = "div", palette="RdBu", direction = 1)

Show normal summaries

set.seed(99)
test_stations <- baseline_stations %>% sample_n(25) %>% pull(Station)

test_clean_summaries <- aus_observation_clean_summary(test_stations, clean_info) %>% unnest(Clean_Yearly) %>% filter(Year_Clean)

test_clean_summaries %>%
    ggplot() + geom_col(aes(x = Year, y = Mean)) + facet_wrap(~ Station) + 
    annotate("rect", 
         xmin = baseline_start, 
         xmax = baseline_end, 
         ymin = -Inf, ymax = Inf, alpha = 0.1)

Evaluate Baselines

test_observations <- aus_clean_observations(test_stations, clean_info)

valuation_info <- aus_valuation_info(
    "RainAnamoly", clean_info, baseline_start = baseline_start, baseline_end = baseline_end, 
    minor_period = aus_monthly_period(),
    input_variables = list(),
    mutators = list(
    ),
    summarisers = list2(
        Monthly_Rainfall = quo(sum(Rainfall, na.rm = TRUE)),
        Sqrt_Monthly_Rainfall = quo(sqrt(Monthly_Rainfall+0.1))
    ),
    output_variable = sym("Sqrt_Monthly_Rainfall"),
    scale = TRUE)

test_valuation <- aus_valuations(test_observations, valuation_info)
test_valuation

test_baselines <- aus_valuation_baselines(test_observations, valuation_info)
test_baselines

test_baselines %>% ggplot() + geom_line(aes(x = Minor_Period, y = Baseline.Mean)) + facet_wrap(~ Station) + ggtitle("Test Baselines means")

test_baselines %>% ggplot() + geom_col(aes(x = Minor_Period, y = Baseline.N)) + facet_wrap(~ Station) + ggtitle("Test Baselines Observations")

plot_baseline_v_time <- function(valuations, station) {
    test_baselines %>% 
        filter(Station == station) %>% 
        ggplot(aes(x = Minor_Period, y = Baseline.Mean)) + 
        geom_line() + 
        geom_errorbar(aes(ymin=Baseline.Mean-Baseline.SD, ymax=Baseline.Mean+Baseline.SD), width=.1) +
        geom_point(aes(x = Minor_Period, y = Baseline.Mean)) +
        ggtitle(station)
}

for (station in test_stations) {
    print(plot_baseline_v_time(test_valuations, station))
}

View anaomolies

test_valuations <- aus_valuations(test_observations, valuation_info, baselines = test_baselines)
test_valuations

test_valuations %>% ggplot() + 
    geom_line(aes(x = Median_Date, y = Rainfall.Anamoly)) + facet_wrap(~ Station) + ggtitle("Test Baselines means")

plot_anomly_v_time <- function(valuations, station) {
    valuations %>% 
        filter(Station == station) %>% 
        ggplot() + 
        geom_line(aes(x = Median_Date, y = Rainfall.Anamoly)) + 
        ggtitle(station)
}

for (station in test_stations) {
    print(plot_anomly_v_time(test_valuations, station))
}

Density plots

test_valuations %>% 
    ggplot() + 
    geom_density(aes(x = Rainfall.Anamoly, color = Station)) + 
    ggtitle(station)

plot_anomly_density <- function(valuations, station) {
    valuations %>% 
        filter(Station == station) %>% 
        ggplot() + 
        geom_density(aes(x = Rainfall.Anamoly)) + 
        ggtitle(station)
}

for (station in test_stations) {
    print(plot_anomly_density(test_valuations, station))
}

Combined plot

test_valuations %>% 
    ggplot() + 
    geom_point(aes(x = Median_Date, y = Rainfall.Anamoly)) 

Full Data

observations <- aus_clean_observations(baseline_stations$Station, clean_info)
observations

baselines <- aus_valuation_baselines(observations, valuation_info)
baselines

valuations <- aus_valuations(observations, valuation_info, baselines = baselines)
valuations

valuations %>%
    filter(Major_Period == 2017, Minor_Period == 1) %>%
    inner_join(baseline_stations, by = c("Zone", "Station")) %>%
    aus_dark_map() + 
    geom_point(aes(x = Longitude, y = Latitude, color = Rainfall.Anamoly)) +
    scale_colour_distiller(type = "div", palette="RdBu", direction = 1)

valuations %>%
    inner_join(baseline_stations, by = c("Zone", "Station")) %>%
    filter(Major_Period > 1990, Longitude > 140, Latitude < -30) %>%
    group_by(Median_Date) %>%
    summarise(Rainfall.Anamoly = mean(Rainfall.Anamoly)) %>%
    ggplot(aes(x = Median_Date, y = Rainfall.Anamoly)) +
    geom_line() +
    geom_smooth()

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'